speizer_vulnerability_testing

library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(readxl)
library(leaflet)


options(
  tigris_class = "sf"
)

Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")

Mapping population of ages 65 and older for the Bay Area

# If not already done, obtain the census data
# acs_vars <-
#   listCensusMetadata(
#     name = "2018/acs/acs5",
#     type = "variables"
#   )
# # Save an .rds version on my computer
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
# saveRDS(acs_vars, file = "censusData2018_acs_acs5.rds")
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")

# obtain the saved census data 
setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
acs_vars = readRDS("censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")

# define Bay Area county names
bay_county_names <-
  c(
    "Alameda",
    "Contra Costa",
    "Marin",
    "Napa",
    "San Francisco",
    "San Mateo",
    "Santa Clara",
    "Solano",
    "Sonoma"
  )

# get Bay Area tracts
bay_tracts <-
  bay_county_names %>% 
  map(function(x){
    tracts("CA",x,progress_bar=F) %>% 
      pull(GEOID)
  }) %>% unlist()

# get sex by age for Bay Area counties, modified from the rBasics.Rmd code
bayArea_sexbyage <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "tract:*", 
    regionin = "state:06",
    vars = "group(B01001)"
  ) %>%
  mutate(
    tractval =
      paste0(state,county,tract)
  ) %>% 
  filter(tractval %in% bay_tracts) %>%
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
  dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
  gather(
    key = "variable",
    value = "estimate",
    -tractval
  ) %>%
  mutate(
    label = acs_vars$label[match(variable,acs_vars$name)]
  ) %>% 
  select(-variable) %>% 
  separate(
    label,
    into = c(NA,NA,"sex","age"),
    sep = "!!"
  )


# Now find the population that is age 65 and up - also modified from rbasics.Rmd
bayArea_elderly <- 
  bayArea_sexbyage %>% 
  filter(!is.na(age)) %>% 
  mutate(
    elderly = 
      ifelse(
        age %in% c(
          "65 and 66 years",
          "67 to 69 years",
          "70 to 74 years",
          "75 to 79 years",
          "80 to 84 years",
          "85 years and over"
        ),
        estimate,
        NA
      )
  ) %>% 
  group_by(tractval) %>% 
  summarize(
    elderly = sum(elderly, na.rm = T),
    total_population = sum(estimate, na.rm = T)
  ) %>% 
  mutate(
    percent_elderly = elderly/total_population
  )


# get Bay Area tracts objects
bay_tracts_obj <- NULL
for (i in 1:length(bay_county_names)) {
  bay_tracts_obj <- bay_tracts_obj %>% rbind(tracts("CA",bay_county_names[i],progress_bar=F))
}


# plot percent elderly
bayArea_elderly %>% 
  left_join(bay_tracts_obj %>% dplyr::select(GEOID), by = c("tractval" = "GEOID")) %>% 
  st_as_sf() %>% 
  filter(!is.na(percent_elderly)) %>% 
  mapview(zcol = "percent_elderly")

Note that the map is very skewed by a few places with very high percent of elderly citizens (one is near the senior living community Rossmoor, for example). Let’s change the legend to make all locations with more than 60% residents 65 and over to show up as grayed out, to be able to see the finer details of everywhere else.

bayArea_elderly %>% 
  left_join(bay_tracts_obj %>% dplyr::select(GEOID), by = c("tractval" = "GEOID")) %>% 
  st_as_sf() %>% 
  filter(!is.na(percent_elderly)) %>% 
  mapview(zcol = "percent_elderly", at = seq(0, 0.6, by = 0.05))

Mapping bed count in the Bay Area

# get data on hospital beds
setwd("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/CHHS")
hospBeds <- read_excel("healthcare_facility_beds.xlsx")
# saveRDS(acs_vars, file = "censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")

# select Bay Area counties
hospBeds <- hospBeds %>% filter(COUNTY_NAME %in% toupper(bay_county_names))

# number of beds at each hospital classified as GENERAL ACUTE CARE HOSPITAL, categorized by facility and type of bed
genAcuteBeds <- hospBeds %>% filter(FAC_FDR == 'GENERAL ACUTE CARE HOSPITAL') %>%
  group_by(FACNAME, COUNTY_NAME, BED_CAPACITY_TYPE) %>% summarize(genAcuteBedsHospAndType = sum(BED_CAPACITY, na.rm = T)) 

# number of beds at each hospital 
genAcuteBedsByHospital <- genAcuteBeds %>% ungroup() %>% group_by(FACNAME, COUNTY_NAME) %>% 
  summarize(totalAcuteBedsPerHosp = sum(genAcuteBedsHospAndType, na.rm = T))

# number of beds in each county
genAcuteBedsCounty <- genAcuteBedsByHospital %>% ungroup() %>% group_by(COUNTY_NAME) %>%
  summarize(totalHosp = n(), totalAcuteHospitalBeds = sum(totalAcuteBedsPerHosp, na.rm = T))

# get geometries for mapview
ca_counties <- counties("CA", cb=F, progress_bar = FALSE)
bay_area_counties_geom <-
  ca_counties %>%
  filter(NAME %in% bay_county_names) %>% 
  mutate(uppercasename = toupper(NAME))

# join
bay_area_counties_geom <- bay_area_counties_geom %>%
  left_join(genAcuteBedsCounty, by = c('uppercasename' = 'COUNTY_NAME')) 

# map total acute care hospital beds by county
bay_area_counties_geom %>% dplyr::select(totalAcuteHospitalBeds) %>% 
  mapview(layer.name = "Total Acute Care Hospital Beds") 

It’s likely more useful to look at the map of beds per capita. Let’s map that as well.

# get the total population in each county
totalPopByCounty <- NULL
for (countyName in bay_county_names) {
  thisCountyTracts <- tracts("CA", countyName, progress_bar=F) %>% 
    mutate(tractval = paste0(STATEFP,COUNTYFP,TRACTCE))
  
  thisCountyPop <- 0
  for (i in 1:nrow(bayArea_elderly)) {
    if (bayArea_elderly$tractval[i] %in% thisCountyTracts$tractval) {
      thisCountyPop <- thisCountyPop + bayArea_elderly$total_population[i]
    }
  }
  totalPopByCounty <- rbind(totalPopByCounty, c(countyName, thisCountyPop))
}

totalPopByCounty <- as.data.frame(totalPopByCounty) 

totalPopByCounty <- totalPopByCounty %>% 
  rename('county' = colnames(totalPopByCounty)[1], 'population' = colnames(totalPopByCounty)[2])

# normalize hospital beds by population, giving in units of beds per 1000 people
genAcuteBedsCountyPer1000 <- bay_area_counties_geom %>% 
  left_join(totalPopByCounty, by = c('NAME' = 'county' )) %>% 
  mutate(bedsPer1000 = 1000*totalAcuteHospitalBeds/as.numeric(levels(population))[population])

# map
genAcuteBedsCountyPer1000 %>% dplyr::select(bedsPer1000) %>% 
  mapview(layer.name = "Acute Care Hospital Beds Per 1000 People", at = seq(1, 3.5, by = 0.5)) %>% 
  addStaticLabels(label = round(genAcuteBedsCountyPer1000$bedsPer1000, digits = 3)) 

These numbers include some beds that may not realistically be available or useful, like intensive care newborn beds. It might be useful to select specifically beds designated for certain uses (e.g. unspecified general acute care, rehabilitation, intermediate care, intensive care, coronary care, and acute respiratory care), assuming these are the ones likely able to be available quickly for COVID-19 patients.